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Abstract 

In an excitable Delaunay triangulation every node takes three states (resting, excited 
and refractory) and updates its state in discrete time depending on a ratio of excited 
neighbours. All nodes update their states in parallel. By varying excitability of nodes 
we produce a range of phenomena, including reflection of excitation wave from edge 
of triangulation, backfire of excitation, branching clusters of excitation and localized 
excitation domains. Our findings contribute to studies of propagating perturbations 
and waves in non-crystalline substrates. 

Keywords: Delaunay triangulation, excitation, waves, localisations, space-time dy- 
namics, pattern formation 



1 Introduction 

Given a finite set of planar points Delaunay triangulation is a planar proximity 
graph which subdivides the space onto triangles with nodes in the given set 
such that the circumcircle of any triangle contains no points of the given 
set other than the triangle's vertices [9]. Delaunay triangulation is a graph- 
theoretic dual of Voronoi diagrams [27j. It represents connectivity of Voronoi 
cells. Voronoi diagram and its dual Delaunay triangulation are widely used 
in studies related to filling a space with connected structural units. Voronoi 
diagram and Delaunay triangulation are used to approximate arrangements 
of discs [H] , sphere packing [20lll()ll21j , to make structural analysis of liquids 
and gases [3j, and protein structure [24], and to model dense gels [28] and 
inter- atomic bonds [T6]. 

We are interested in studying excitable Delaunay triangulation because they 
may provide a good alternative to existing approaches of modelling unstruc- 
tured unconventional computers [25j . Experimental research in novel and emerg- 
ing computing paradigms and materials shows a great progress in designing 
laboratory prototypes of spatially extended computing devices. In these de- 
vices computation is implemented by excitation waves and localisations in 
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reaction-diffusion chemical media [2], geometrically constrained and compart- 
mentalized excitable substrates p!2l[T3l[8l[T9] . organic molecular assemblies J5|, 
and gas-discharge systems [4j . These unconventional computing substrate can 
be formally represented by Delaunay triangulations with excitable nodes. Thus 
it is important to uncover most common types of excitation dynamics on the 
Delaunay diagrams. 

Despite being a ubiquitous graph representation of wide range of natural phe- 
nomena the Delaunay triangulation was not studied from automaton point 
of view. Will excitable Delaunay triangulation behave as a conventional ex- 
citable cellular automata or there will be some unusual phenomena? We an- 
swer the question by slightly modifying classical Greenberg-Hasting model [15] 
and considering not only a threshold of excitation but also a ratio of excited 
neighbours as an essential factor of nodes' activation. 

The paper is structured as follows. We introduce automata on triangulations 
and excitation rules in Sect. [2| In Sect. |3] we discuss structural properties of 
automata triangulations. Sections [4] and [5] present classification of space-time 
dynamics of excitation for absolute (based on a number of excited neighbours) 
and relative (based on a ratio of excited neighbours) rules of excitation. Results 
are discussed in Sect. [HI 



2 Delaunay automata and excitations 

Given a planar finite set V the Delaunay triangulation [9] ^^(V) = (V, E) is 
a graph subdividing the space onto triangles with vertices in V and edges in 
E where the circumcircle of any triangle contains no points of V other than 
its vertices. Neighbours of a node G V are nodes from V connected with v 
by edges from E. 

The set V is constructed as follows. We take a disc-container of radius 480 
and fill it with up to 15,000 disc-nodes. We assume that each disc-node has 
radius 2.5, thus a minimal distance between any two nodes is 5. The Voronoi 
diagram, and its dual triangulation, are appropriate representations of such 
identical sphere packing on 2D surface, where planar points of V represent 
centres of the spheres. 

We define density as a ratio of areas occupied by discs to area of the disc 
container. Examples of triangulations for densities cp = 0.0027,0.027,0.136 
and 0.407, corresponding to number of nodes packed 100, 1000, 5000, and 
15000, are shown in Fig. [T] 

The Delaunay automaton is defined in the following way. A node G V is a 
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finite state maciiine. Every node updates its state in discrete time depending 
on states of its neiglibours. All nodes update their states simultaneously. Nodes 
can have different number of neighbours therefore we better use totalistic node- 
state update function, where a node updates its state depending on just the 
numbers of different node-states in its neghbourhood. 



Here we are concerned only with modelling excitation on Delaunay triangula- 
tions. Thus we assign three states — resting (o), excited (+) and refractory 
(— ) — to nodes of V. We assume that a resting node excites depending on 
a number of excited neighbours. If a node is excited at time t the node takes 
refractory state at time step t + 1, independently on states of its neighbours. 
Transition from refractory to resting state is also unconditional. 



Let u{v) = {u : (vu) G E} be node '^'s neighbourhood, a state of node 
V at time step t, a^{v) a number of excited neighbours of v at step t, and d{v) 
be a degree, or a number of neighbours |^^(^)|, of node v. Then the node-state 
transition functions can be defined as follows. 



• Absolute excitability: 



iia\v) G S 
if = + 
o, ii — — 



(1) 



Node V excites if a^{v) G S, where S is a set of natural numbers. For 
example, S = {2,4} means a resting node excites if it has two or four 
excited neighbours. The rule includes threshold excitation a^{v) > 9^ 9 is Si 
natural number. For example, 9 = 2 means a resting node excites if it has 
at least 2 excited neighbours, i.e. S = {2, 3, 4, • • • }. 
• Relative excitability: 



+, if ^ >6 

' d(v) 
if = + 

o, if = — 



(2) 



Node V excites if p^{v) = > e, where < 6 < 1. This condition bring 
more 'fairness' in excitation process. Some nodes can have less neighbours 
than other nodes, thus measuring excitation of neighbourhood just by num- 
ber of excited neighbour would not be 'fair'. It is feasible to calculate a ratio 
p^{v) of excited neighbours a^{v) to a total number of neighbours d{v). 
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3 Structural properties of Delaunay automata 



Propagation of excitation in Delaunay automata depends on structural prop- 
erties of their graphs. We found that in case of sparse distribution of disc-nodes 
(Fig. [T^), (j) = 0.0027, majority of nodes have five then six and four neigh- 
bours each (Fig. [2^). With increase of the density (Fig. [l|3-d) maximum of 
degree distribution is shifted to six neighbours per node (with probability 
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Fig. 3. Dynamics of occupation a of triangulations from single-node perturbation; 
a is a ratio of occupied nodes at time step t to a total number of nodes in the 
triangulation. Graph presents dynamics of occupation for densities (j) — 0.136 and 
= 0.407, and occupation thresholds 77 = 1 and 77 = 2. 

p(6) = 0.45), five nodes {p{h) = 0.27) and seven nodes {p{7) = 0.22) (Fig.|2^). 
The deviation of the distribution significantly decreases with increase of the 
density (Fig. 

The degree distributions found experimentally conform to classical results on 
packing of equal spheres which is between 5.5 for sparse packing and 6.4 for 
a close yet random packing [6], and mean degree of 6 for randomly packed 
sphere with density 0.64 fl4j. 

Let every node of a triangulation takes just two states: '0' (unoccupied) and 
'r (occupied). A resting node becomes occupied if it has at least rj neighbours 
in state T. Initially just one (for = 1) or two (for r] = 2) nodes are assigned 
state 'r. Dynamics of occupation, measured in a ratio a of nodes occupied by 
time step t, is shown in Fig. (S) 

A speed of propagation of state T measured in a ratio of nodes occupied 
in step t increases with decrease of occupation threshold ry. If we measure 
speed in discrete time steps, we will see that it also decreases with increase 
of the density (j). However, the distance covered in any period of time remains 
comparable between triangulations with different number of nodes. 

Finding 1 Threshold rj of node occupation determines a shape of propagating 
occupation fronts. 

For T] = 2 occupation wave-front is convex, shaping into planar by the end of 
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Fig. 4. Time lapsed contours of activation in triangular lattices with density 
(j) — 0.407 and activation threshold (a) 77 = 1 and (b) 77 = 2. Propagating wave 
front is converted to a contour every 10th step of simulation. Initially east most 
nodes are activated. 

propagation (Fig. ^p)- In case of lower threshold of occupation, = 1, prop- 
agation wave near the edges of triangulation moves quicker then inside the 
triangulation core (Fig. [4^). Thus wave front becomes concave. The part of 
wave front travelling along edges of triangulation reaches the side of triangula- 
tion opposite to the initial perturbation side quicker then the front propagating 
inside the triangulation. Thus wave front becomes closed (Fig. [4^). The do- 
main surrounded by a nested group of target waves, in the western part of 
the triangulation in Fig. [4^, is the place where the occupation waves traveling 
from east collapse. 

Finding 2 Perturbation propagates faster along edges of triangulation when 
threshold of node perturbation is low. 

An edge of triangulation is a set of nodes lying on segments of the convex hull 
of V. Spatial structures of node neigbourhoods at the edge of triangulation 
may be responsible for the particulars of propagations described above. 

Typically neighbours of each node are distributed more or less equally around 
the node while nodes belonging to the edge of triangulation have their neigh- 
bours located towards inside of the triangulation or on the edge (Fig. [5^). We 
can integrate spatial distribution of neighbours of node as a vector p = Tfg 
from the node v to geometric centre g of the node v^s neighbourhood. Edge 
nodes of triangulation have usually longer vectors p, see Fig. [5)3. This is why 
a perturbation propagates faster on or near edges of triangulation for a low 
threshold of occupation (Fig. Ilk). Increase of occupation threshold is disad- 




Fig. 5. Geometry of neighbourhoods: (a) a fragment of triangulation, density 
(j) — 0.136, with visible neighbourhoods of edge nodes, edge of triangulation is 
on the right, (b) vectors from nodes to their geometric centres, only vectors which 
length exceeding eight units are shown. Density is = 0.136. 

vantageous for edge nodes which is reflected in changed shape of the growing 
pattern (Fig. ^p)- 



4 Absolute excitability 



If a resting node v excites when it has at least one excited neighbour, rule ([T]) 
a^{y) > 1, 'classical' excitation waves are observed (Fig.[6^b). By 'classical' we 
mean that excitation waves annihilate when they reach edges of triangulation, 
two waves merge when they collide one with another, a single-site perturbation 
initiates a circularly propagating excitation wave, and refractory state is a 
necessary component in seeds of target- wave generators. 

The excitable triangulation behaves less conventionally when nodes are selec- 
tive in their excitability. Consider the rule ([T]) S = {1}: a resting node excites 
if it has exactly one excited neighbour. Single site excitation leads to forma- 
tion of circular waves. The waves do merge when collide with each other and 
also they may form generators of target waves at the sites of their collision 
(Fig. 15). 
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Fig. 6. Examples of excitation dynamics for rule ([l]) g^{v) > 1 (a)-(c) and rule ([l]) 
a^(v) = 1 (d): (ab) waves generated by singular excitations merge when collide, 

(b) random initial configuration (at the beginning of simulation a node got excited 
state with probability 0.1, refractory state with probability 0.1, and resting state 
with probability 0.8) develops into a configuration of several target-wave generators; 

(c) configuration of excitation developed in 300 time steps from initial configura- 
tion where all but one nodes are resting. Density of disc-nodes in triangulation is 
(j) = 0.407. 

Finding 3 Let a resting node of Delaunay triangulation excite if exactly one 
neighbour is excited^ then a generator of target-waves can be produced by a 
single-site excitation. 

For a^(v) > 2 no excitation persists. However when resting node excites if 
exactly two or three of its neighbours are excited we obtain a quasi-chaotic 
excitation dynamics (Fig. |6]l). This happens when nodes follow state update 
rules g for S = 2 and S = {2, 3}. 



1 . 



9 



(a) t = 35 (b) t = 52 




(c) t = 65 (d) t = 148 

Fig. 7. Example of wave generators formed due to backfiring of excitation from 
a single propagating wave front. Rule a^(v) > 1 with delayed (by 10 time steps) 
recovery from refractory state. Density of disc-nodes in triangulation is = 0.407. 

A delay in recovery from refractory states modifies space-time dynamics of 
excitation. Assume that if a node took refractory state then the node stays 
in the refractory state for S time steps. In automata governed by rule ([T]) 
with excitability a^{v) > 1 initial random disturbances — nuclei of excited 
and refractory states — will not lead to formation of wave generators (as in 
rule a^{v) > 1 without node-state recovery delay). However a phenomenon of 
excitation backfiring is still observed. Travelling excitation wave can backfire 
with localized excitations. These localized excitations pass through the wave's 
refractory tail and initiate generation of new wave fronts (Fig. [7]). 

Finding 4 Delaunay excitable automata governed by rules of absolute ex- 
citability exhibit the following phenomena: 
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(a) t = 20 



(b) t = 30 



(c) t = 41 




\ 

(d) t = 67 (e) t = 100 (f) t = 132 



Fig. 8. Dynamics of excitable triangulation for rule ([2]) e = 0.091. Initially just a 
single node is excited. Density of disc-nodes in triangulation is = 0.407. 

• threshold activation rules causes formation of classical excitation wave fronts, 

• rules relying on exact number of excited neighbours show formation of target 
wave generators during collision between ordinary circular waves, 

• when recovery of a node from refractory state to resting state is delayed 
propagating wave fronts backfire with localized excitations, which cause for- 
mation of target-wave generators. 



5 Relative excitation 

We found that we can initiate a persistent excitation patterns for < e < 0.25. 
For e > 0.25 any initially invoked excitation quickly ceases. Automata excited 
by rule ^ with e = exhibit persistent global oscillations because every node 
autonomously follows the cycle o ^ + ^ — ^ • • • . For excitability range 
< 6 < 0.09 Delaunay automata exhibit classical excitation waves. Single ex- 
citation generates single circular wave, clusters of excited and refractory states 
may become generators of target wave. When a propagating wave front hits 
edge of triangulation the wave disappears. Two colliding wave-fronts merge or 
annihilate. 
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(a) t = 18 



(b) t = 28 



(c) t = 37 




(d) t = 45 (e) t = 53 (f ) t = 83 



Fig. 9. Example of backfiring of excitation wave fronts which leads to formation of 
generators of target waves in excitable triangulation, density (j) = 0.407, for rule ^ 
with e 0.11. 

When threshold of excitation increases to e = 0.09 edges of the triangulation 
may become reflective for excitation waves. If an excitation wave front hits an 
edge of the triangulation the front reverses it velocity vector and again propa- 
gates inside the triangulation. An example is shown in Fig. [8| Initially graph is 
in resting state. We excited one node. A circular wave of excitation propagates 
outwards the initial perturbation (Fig. [8^). When the wave front reaches edge 
of the triangulation it starts to disappear and almost annihilate. However in 
some part of triangulation edge domains of centrifugal wave-fragments evoke 
centripetal wave-fragments (Fig.|8}3c). Fragments of the centripetal wave closer 
to edge of triangulation propagate quicker then those inside the triangulation. 
Therefore the wave encircles the triangulation (Fig. [sjle) and collapses inside 
the triangulation (Fig. [sjf). 

The reflection of waves can be observed for values 0.09 < e < 0.11. Further 
increase of excitation threshold brings up the phenomenon of excitation back- 
firing again. 

Triangulations, which nodes are excited if ratio e is at least 0.11, exhibit back- 
firing — generators of target waves are formed behind the wave front initiated 
by a single excitation (Fig. [9]). Due to inhomogeneous structures of node neigh- 
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Fig. 10. Example of generators formed close to a site of original perturbation. Au- 
tomaton is excited by rule ^ with e = 0.150. Density of disc- nodes in triangulation 
is = 0.407. 



bourhoods the wave front (Fig. [9^) breaks up (Fig. [9p). A gap is formed and 
some part of the front folds backward. A spiral wave or waves are formed. 
They give rise to a succession of target waves (Fig. [oj^de). Eventually origi- 
nal wave front disappears at the edge of triangulation but the triangulation 
remains filled with sources of target waves (Fig. |9]f). 

With 6 increasing from 0.11 to 0.17 a number of wave generators forming 
behind propagating wave front increases considerably. Thus, in networks ex- 
cited by rule ^ with e > 0.14 a generator is formed almost immediately after 
single-excitation wave starts its propagation. The earlier in time a generator is 
born the more likely the generator will dominate the triangulation. Therefore 
in many cases we can observe just few generators close to the site of original 



excitation (Fig. 10). 



With e exceeding 0.17 the automaton approaches regime of sub-excitability. 
Single site excitation no longer leads to formation of a circular wave. However 
we observe travelling and stationary localized excitations, distant analogous 
to wave- fragments in sub-excitable Belousov-Zhabotinsky medium [7J. Initial 
local perturbations lead to propagating localized excitations. The excitations 
later can form a localized, and not changing it is outer shape, domain of activ- 



ity, see example in Fig. [11} or a slowly growing domain, which may eventually 
occupy the whole triangulation (Fig.[l2]). Usually, two or three wave-fragments 
are formed near the site of initial activation of resting triangulation. These 
wave-fragments travel outward the perturbation loci. Some time after their 
initiation the wave-fragment backfires few excited micro-localizations which 
may form generators of wave-fragments. Thus the domain of excitation is 
born. Sometimes generators of wave-fragments emerge due to folding of the 
'wings' of a wave-segment backwards. Folded parts of wave-fragment interact 
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(a) t = 10 




(b) t = 16 




(c) t = 40 (d) t = 142 

Fig. 11. Example of excitation dynamics for e = 0.17. A local perturbation leads 
to localised excitations and slowly growing domains. The domain shown in (d) has 
stationary boundaries however configuration of excited and refractory states inside 
the domain is changing. Density of disc-nodes in triangulation is = 0.407. 

with each other and thus they produce the generator. 



The growing domains guided by mobile localisations are typical for excitability 
values 0.17 < 6 < 0.2. When e exceed 0.2 the domains seize propagating. 
The excitation can still persist but in a minuscule quantity. A random initial 
configuration, where every node gets one of three states equiprobably, is either 
transformed to a totally resting state or to a resting configuration with one 
or two tiny oscillators. Examples of most common oscillators are shown in 



Fig. 13 



What are degrees of nodes occupied by non-resting states of the oscillators? 
We stimulated triangulations with random configurations of excited and re- 
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(d) t = 53 (e) t = 119 (f) t = 278 



Fig. 12. In triangulation, density of disc-nodes is = 0.407, excited by rule ^ with 
e = 0.17 a local perturbation leads to localised excitations and slowly propagating 
domains, which may occupy significant number of nodes. 




(c) 



Fig. 13. Examples of oscillators: (a) and (b) two oscillators for excitability e = 0.2, 
(c) minimal oscillator, e = 0.21. Each sub-figure presents three consecutive states 
of an oscillator. Density of disc-nodes in triangulation is = 0.407. 
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0.112 



0.25 



Fig. 14. Visualisation of basic types of excitation activity parameterised by relative 
excitability threshold e. 

fractory states (each state is assigned with probability |), waited till all activ- 
ity but minimal oscillator ceases and recorded degrees of nodes occupied by 
the oscillator states. We found that nodes occupied by excited or refractory 
states in the minimal oscillator (Fig. have 4, 6, 7 or 8 neighbours. In 
larger oscillators nodes have degrees 4,5,6,7,8 and 4,5,6,7,8,9. 

Proposition 1 Diversity of node degrees is a necessary requirement for exis- 
tence of an oscillator in a sub-excitable triangulation. 

The oscillators survive for 0.2 < e < 0.25. For e > 0.25 no excitation persist 
at all. 

Finding 5 Delaunay excitable automata governed by rules of relative excitabil- 
ity with dimensionless threshold of excitation e exhibit the following phenom- 



ena: 




e 


phenomenon 





global oscillations 


< e < 0.09 


classical excitation waves 


0.09 < e < 0.11 


waves are reflected from edges 


0.11 < e < 0.17 


waves backfire with excitation 


0.17 < 6 < 0.2 


localized wave-fragments lead growing tips of branching domains 


0.2 < e < 0.25 


only tiny oscillating domains are present 


0.25 < e 


no excitation persists 



The finding is illustrated in Fig. [Mj Is there any structural meaning of e- 
boundaries between different classes of excitable triangulations? Value e = 
0.09 corresponds to a situation when a node with 11 neighbours has at least 
one excited neighbour. Nodes with degree 11 are rare species in Delaunay 
triangulations (Fig. [2]). However when such nodes are deprived, i.e. for e > 
0.09, from a chance to be excited at all by a minimal perturbation waves of 
excitation start to be reflected by edges of triangulation. 
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Excitation wave- fronts backfire when we make nodes with degree 9 non-excitable 
by raising excitabihty to over 0.11. This is because e = 0.11 describes a situ- 
ation of a node with 9 neighbours which has just one neighbour excited. 

Locahzations emerge in triangulations when e > 0.17. In principle this value 
of e corresponds to two situations: a node with six neighbours has one ex- 
cited neighbour and a node with 12 neighbours has two excited neighbours. 
The latter situation is not important because nodes with 12 neighbours are 
extremely rare. Thus, we can propose that excitation becomes localized when 
six-neighbour nodes are activated by at least two neighbours. Two excited 
neighbours is a minimal requirement for existence of localizations in orthogo- 
nal automaton networks with eight-cell neighbourhoods In a regular net- 
work, or cellular automaton, a travelling localization preserves its shape and 
velocity vector as long as necessary. 

Delaunay triangulations are irregular therefore localized wave-fragments do 
not survive for a long time. These travelling localizations frequently change 
directions of their motion. Soon after its birth a localization either perishes or 



expands, backfires and forms a slowlty growing cluster of activity (Fig. 12). 
Any excitation ceases to sustain when e exceeds 0.25, which corresponds to 
one of four, two of eight, and three of twelve excited neighbours. 

Does density of disc-nodes affect structure of behavioural space? 

For relative excitability < e < 0.09 decrease of density increases chances 
of wave- generators to be formed from a simple propagating wave- front. This 
possibly happens due to increased inhomogeneity of node degrees and in- 
creasing role of each particular node in excitation propagation. In case of 
0.09 < e < 0.11 decrease of cj) increases opportunities for wave-generators to 
form when a wave collides to edge of triangulation. 

For excitability 0.11 < e < 0.17 time of backfiring and formation of genera- 
tors behind wave-fronts, is proportional to density 0, higher the (j) the longer 
it takes for a singular wave- front to backfire. And finally, chances that an 
excitation domain growing from a single perturbation, in triangulations with 
excitability 0.17 < e < 0.2, occupies the whole triangulation are inversely 
proportional to density of disc-nodes in the triangulation. 



6 Discussion 



We defined excitable automata on Delaunay triangulation and demonstrated 
how to control a space-time dynamics of excitation on the triangulation using 
absolute and relative excitability thresholds. We uncovered several interesting 
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phenomena ranging from reflection of excitation waves by edge of triangulation 
to branching domains of activity guided by traveUing locahsed excitations. 

Growing domains of locahzed excitation may be analogous to patterns of exci- 
tation in Belousov-Zhabotinsky reaction in micro-emulsion [18]. Our flndings 
on reflection of waves by edge of triangulation support previously published 
results on active waves reflection in spatially inhomogeneous non-linear me- 
dia and reflection of reaction-diffusion waves by no-flux boundary due to 
consumption of reactant between a wave and a boundary [23j . Further studies 
of wave speed up along edge of triangulation could be enhanced by additional 
techniques of characterising topological and dynamical properties of complex 
networks, e.g. a diversity entropy proposed in [26j. 

We believe our flndings will contribute towards designs of novel computing sub- 
strates in non- crystalline structure. Also, automaton interpretation of activity 
dynamics on Delaunay triangulation can make a viable model of automaton- 
network approaches to design of nano-computing devices [17j. 
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